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Recently, a lot of effort has been put into describing the thermalization of the quark-gluon plasma 
using the gauge/gravity duality. In this context we here present a full numerical solution of the 
early far-from-equilibrium formation of the plasma, which is expanding radially in the transverse 
plane and is boost invariant along the collision axis. This can model the early stage of a head-on 
relativistic heavy ion collision. The resulting momentum distribution quickly reaches local equi- 
librium, after which they can be evolved using ordinary hydrodynamics. We comment on general 
implications for these hydrodynamic simulations, both for central and non-central collisions, and 
including fluctuations in the initial state. 



1. Introduction. Describing the very early stage of a 
relativistic heavy ion collision has remained challenging, 
in particular when the state is far-from-equilbrium. In 
this stage QCD is partially strongly coupled and impos- 
sible to solve in general. However, a lot of progress has 
been made by employing a dual gravitational description 
(see pQ and its references), the two seminal examples be- 
ing the fast thermalization [21 [3] and the very small vis- 
cosity in the consequent hydrodynamic regime [3]. 

These dual descriptions have been studied for near- 
equilibrium physics 4-6 and also (numerically) far-from- 
equilibrium [51 171-H^]. Yet, these studies have always as- 
sumed homogeneity in the transverse plane, which in par- 
ticular makes it impossible to study radial flow. Recently, 
Bantilan, Gubser and Pretorius [T3] have performed a 
simulation having both longitudinal and radial expan- 
sion. However, in their set-up it was not possible to have 
boost-invariance in the longitudinal direction and, more 
importantly, the evolution was always in the hydrody- 
namic regime. In this Letter we will present a numeri- 
cal simulation having radial flow in the transverse plane, 
boost-invariance in the longitudinal direction, and being 
far-from-equilibrium initially. 

Essentially, we follow the numerical scheme worked out 
in [5], but we assume boost-invariance in the longitudi- 
nal direction and allow for non-trivial radial dynamics in 
the transverse plane. This keeps the gravitational prob- 
lem 2+1 dimensional, which can be solved using pseudo- 
spectral evolution. As initial conditions we present two 
simple models; the first starts with a blob of energy with 
a diameter of approximately 14 fm in vacuum, whereas 
the second has a blob of about 1 fm in a bath of half 
the peak energy density. These initial states can model 
the overall thermalization of a central collision and the 
evolution of an initial fluctuation in such a collision (fluc- 
tuations recently became a topic of much interest, see 
for instance [H]). For the bulk metric we started with 
vacuum AdS, but adapted the near-boundary coefficients 
for the energy density and the pressures according to the 
Glauber model. 

The results of our simulations are both intuitive and 
encouraging. Firstly, we find that our geometry thermal- 



izes very quickly, confirming previous studies. Secondly, 
our radial velocity profile at the end of our evolution is 
similar to typical initial conditions used for a hydrody- 
namical evolution. So these radial velocities serve as a 
confirmation that current hydrodynamic simulations do 
not have to be modified dramatically, but they also pro- 
vide an improvement for hydrodynamic evolutions. 
2. Holographic model. As our coordinates in the field 
theory it is natural to use proper time r and rapidity y, 
defined by t — rcoshy and x\\ = rsinhy, and angular 
coordinates p and 9 in the transverse plane. The assump- 
tions of boost-invariance and rotational symmetry then 
imply that all functions are independent of y and 8. In 
these coordinates the flat metric of the field theory reads 



ds% = -dr 2 + dp 2 + p 2 d0 2 + r 2 dy 2 



(1) 



Given these symmetries and using generalized 
Eddington-Finkelstein coordinates, we can write 
the dual bulk metric as 

ds 2 = -Adr 2 + Y?{e- B - C dy 2 + e B dp 2 + e c d0 2 ) 

+2drdr + 2Fdpdr, (2) 

where A, B, C, S and F are all functions of r, p and the 
bulk radial coordinate r. Solving Einstein's equations 
(with cosmological constant A = —6) order by order in 
r, demanding that ds 2 \ r=QO = r 2 ds 2 B , gives us the near- 
boundary expansion of the bulk metric: 
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where in these expressions we fixed a residual gauge free- 
dom r — > r + £(t, p) by demanding d r A\ r=OQ = 2r. 
The normalizable modes of the metric, a 4 , 64, C4 and 
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fi, depend on the full bulk geometry and cannot be de- 
termined from a near-boundary expansion. Using the 
gauge/gravity duality we can now determine the stress 
tensor of the dual field theory [T5] , which has five inde- 
pendent non-zero components: 
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all functions of r and p, where we put the number of 
colors N c = 3. The conservation of the stress tensor 
implies that 
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Our model basically contains two scales: the initial en- 
ergy density and the characteristic scale in the radial 
direction. We can, however, make use of the scale in- 
variance of the field theory to rescale our coordinates 
such that at t = 0.6 fm the energy density at the ori- 
gin equals £0 = 187GeV/fm 3 (3TJ. We choose this com- 
bination to reproduce the final multiplicities of central 
heavy- ion collisions at LHC [T6] . For the radial profile 
we then consider two types of initial conditions, speci- 
fied at some small time Tj n ss 0.12 fm [52]. The first is 
a model for a head-on collision, where the shape of the 
energy density is provided by the Glauber model, having 
an approximate radius of 6.5 fm. The second energy den- 
sity profile models one fluctuation in the initial state of 
such a collision. We take a Gaussian of width 0.5 fm for 
this profile (see figure [T]). For both initial conditions we 
assume that initially there is no radial momentum, such 
that / 4 (T m , p) = 0. 

Importantly, we must also specify the metric functions 
B(r, Ti n , p) and C(r, T% n , p) on a full time-slice of the 
bulk AdS geometry. These two functions, together with 
04, fi and the Einstein equations, specify the complete 
metric and its time derivative on a time-slice [8] . In prin- 
ciple, these functions should follow from a model describ- 
ing the very first weakly coupled stage after the collision, 
such as the Glauber model or the Color Glass Conden- 
sate. However, these models themselves contain signifi- 
cant uncertainties and, more importantly, it is not clear 
how to map them to this gravitational setting. Therefore 
we made a simple choice, where B and C are the same 
functions as in vacuum AdS, but with modified 64 and C4, 
such that the longitudinal pressure p y vanishes initially. 




FIG. 1: The initial energy density profiles at Ti n = 0.12 fm 
as a function of the distance to the origin. The blue curve 
models a central heavy-ion collision; the red curve models a 
fluctuation in such a collision. 



Having specified the initial and boundary conditions 
we can solve Einstein's equations numerically [23], using 
essentially the same scheme as in [8 . One small difference 
is the required boundary conditions in the p direction, 
which in this case means smoothness at the origin and at 
infinity. As in [5], we added a small (3%) regulator energy 
density and checked that our results do not depend on 
this regulator. 

3. Results. After determining the stress tensor one 
can extract the radial velocity, defined by the boost af- 
ter which there is no momentum flow. Figure [3] shows 
this velocity times the energy density, which gives a good 
measure of the momentum flow. The radial velocity, to- 
gether with the stress tensor in the local rest frame, can 
be used to compute the stress tensor according to hy- 
drodynamics. Although initially there will not be local 
equilibrium, at late times a hydrodynamic expansion is 
expected to be valid. It is therefore interesting to com- 
pare the actual pressures with the pressures which follow 
from a hydrodynamic expansion [171 118] . 

In figure [4] we plot the difference of p p and the cor- 
responding first order hydrodynamic prediction of our 
model of a nucleus. The stress tensor is excellently de- 
scribed by hydrodynamics as soon as r = 0.35 fm. At 
the border of our nucleus this is slightly subtler, since 
the stress tensor is rather small there, and it becomes 
comparable to our regulator energy density. We there- 
fore cannot say too much about this, but the agreement 
with hydrodynamics is also there encouraging. We note 
that in previous studies jTDl [TTJ somewhat larger ther- 
malization times (with respect to the local temperature) 
were found, so we expect more exotic initial conditions 
in our bulk AdS to give somewhat later thermalization. 

In figure^ we plot the radial acceleration of our model 
of a fluctuation. We notice the acceleration already de- 
creases considerably during our simulation, in contrast 
with the model for the nucleus. Also, the acceleration 
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FIG. 2: (a) The radial acceleration of our nucleus model. The acceleration decreases after some time, which is mainly a 
consequence of the decrease in radial pressure, due to the isotropization. Thereafter the acceleration is quite steady and mainly 
localized near the boundary of the nucleus, (b) The radial acceleration of our fluctuation model. Since the bump of energy is 
much smaller one can clearly see the spreading out and the decrease in acceleration. As will also be clear from figure [5j this 
model reaches a lower radial speed than the model for the nucleus. 




FIG. 3: The radial velocity times the energy density as a 
function of proper time r and distance to the origin p for 
our model of a nucleus. Note that at late times the increas- 
ing velocity is almost exactly compensated by the decreasing 
energy density (which is due to the longitudinal expansion). 
The slope at the origin at the end of our simulation equals 
0.66 GeV/fm 4 . 




FIG. 4: The difference between the full non-equilibrium p p 
and the pressure given by first order hydrodynamics. Al- 
though hydrodynamics applies very quickly, the viscous con- 
tribution is still large (shown by a red lines). The relatively 
high values for p > 7 fm are a consequence of the very small 
energy density. For the model of a fluctuation the graph is 
similar, with equally quick thermalization. 



increases rapidly near the origin, whereas for the nucleus 
it is rather narrowly peaked near the boundary of the 
nucleus. This means that fluctuations are expected to 
spread out rather quickly. Perhaps surprisingly, also the 
stress tensor for the fluctuation is governed by hydrody- 
namics within 0.35 fm. 

4. Discussion. The main motivation for this study is 
to provide a description of the far-from-equilibrium stage 
of heavy-ion collisions, including non-trivial dynamics in 
the transverse plane. While we kept rotational symme- 
try in the transverse plane, we believe our study can be 
used more generally. One reason for this is an old re- 
sult in asymptotically fiat space [TU], recently studied 
in asymptotically AdS [11], that during black hole for- 
mation gravity can be well approximated by linearizing 
around the final state. We therefore believe that an ini- 
tial energy profile with many fluctuations could be well 
approximated by superposing the result of our fluctua- 



tion presented above. 

Also, it should be possible to use our results for non- 
central collisions. Firstly, one should notice that the ther- 
malization time is much shorter than the size of a nucleus, 
which makes a local analysis valid. Locally the energy 
density and the absolute value of the gradient of the en- 
ergy density are then the most important scales, which 
will determine the (approximate) final velocity at that 
point, at the moment of thermalization. Figure [5] plots 
this relation and one can see that the plots are qualita- 
tively and also quantitatively similar, which confirms the 
argument of this paragraph. This suggests the following 
formula for velocities in the transverse plane: 

Mr) « ^P^eo. (6) 

3e 

This formula would also be the (naive) guess if perfect 
fluid hydrodynamics would apply from the start. Al- 
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FIG. 5: Here we plot the velocity at time r = 0.6 fm (solid 
lines) and r = 0.35 fm (dashed lines), as a function of the 
initial energy density gradient, divided by the energy den- 
sity, for our two models. The two models behave surprisingly 
similar, considering that our two models describe very differ- 
ent dynamics. The red curve lies lower, which is mainly a 
consequence of the decreasing acceleration displayed in fig- 
ure [2|j. We suggest that in a non-symmetric initial state one 
may model the radial velocities as local functions of the initial 
energy density by formula [6] shown as black lines. 

though this is definitely not true, in our cases viscous 
corrections compensate non-hydrodynamic behavior in 
such a way that the formula holds remarkably well. At 
t = 0.35 this is true for both models, whereas after that 
the velocity of the fluctuation is lower, due to the de- 
creasing acceleration (displayed in figure [2]d). 

Apart from the well-known fact that the gauge/gravity 
duality in our setting describes a supersymmetric theory 
at infinite coupling and a large number of colors, there 
is another reason why one should take care applying our 
results directly to experimental settings. Our initial bulk 
metric should, in principle, follow from the way the ex- 
perimental state was created. It is not clear how to do 
this and we therefore adopted a simple model, where the 
metric is basically close to vacuum. In future work we 
plan to study the dependence of the outcomes on this 
initial state, but preliminary findings suggest that the 
results do not depend strongly on the initial bulk metric. 

Even though our model is much simpler than a real 
heavy-ion collision in QCD, we carried out a far-from- 
equilibrium calculation at strong coupling with non- 
trivial dynamics in the transverse plane. These results 
can therefore be very useful as initial conditions for the 
hydrodynamic modeling of heavy-ion collisions. Espe- 
cially the radial velocity at this initial time was basically 
unknown and is hence usually taken to be zero. Our 
radial velocities would be a more natural start for a hy- 
drodynamic simulation. Although the final effect on ex- 
perimental observables would be rather moderate at the 
moment, it will become increasingly important as the ex- 
perimental data improve. 
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